---

Exploring spatial data

Let me get the files we created in Python:

rm(list = ls())
library(sf)

## Warning: package 'sf' was built under R version 4.4.1

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

mainLink='https://github.com/DACSS-CSSmeths/Spatial-Exploring/raw/refs/heads/main/'
mapLink=paste0(mainLink,'maps/countriesCIA.gpkg')


countriesCIA=read_sf(mapLink, layer='cia')
worldBorders=read_sf(mapLink, layer='border')

We have these interesting variables:

names(countriesCIA)

##  [1] "COUNTRY"              "name"                 "region"              
##  [4] "obesityAdults_rate"   "TobaccoUse_perc"      "Alcohol_LitersPerCap"
##  [7] "tobacco_code"         "tobacco_levels"       "tobacco_custom"      
## [10] "alcohol_code"         "alcohol_levels"       "geom"

A Choropleth for CONTINUOUS values

Here you find the last version for this case:

library(ggplot2)

base=ggplot(data = worldBorders)+geom_sf(fill='grey',color=NA) + theme_linedraw()
base + geom_sf(data=countriesCIA,
               aes(fill=TobaccoUse_perc), #variable for coloring geometry
               color=NA) + # no borders
    labs(fill="Tobacco use\n(grey = missing)") +
    scale_fill_viridis_c(direction = 1,option='A') # this is MAGMA!

A Choropleth for DISCRETIZED values

Let’s plot a diverging palette (choose from here) using the discrete version:

base + geom_sf(data=countriesCIA,
               aes(fill=tobacco_levels),color=NA) + 
               labs(fill="Level",
                    title='Tobaco use')+
                    scale_fill_brewer(palette = "PiYG",direction = -1) 

Creating layers from the discrete values

customCols=c("green", "orange")
base + geom_sf(data=countriesCIA[countriesCIA$tobacco_code%in%c(0,2),],
               aes(fill=tobacco_levels),color=NA) + 
               labs(fill="Level",
                    title='Tobaco use')+
                    scale_fill_manual(values = customCols)+
    facet_grid(~tobacco_levels) + guides(fill="none")

Beyond choroplething

interestingCountries=countriesCIA[(countriesCIA$tobacco_code==4) & (countriesCIA$alcohol_code==4),]
interestingCountries

## Simple feature collection with 4 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 13.50479 ymin: 34.64027 xmax: 34.58604 ymax: 58.08326
## Geodetic CRS:  WGS 84
## # A tibble: 4 × 12
##   COUNTRY  name   region obesityAdults_rate TobaccoUse_perc Alcohol_LitersPerCap
##   <chr>    <chr>  <chr>               <dbl>           <dbl>                <dbl>
## 1 Bulgaria Bulga… Europe               25              39                  11.2 
## 2 Cyprus   Cyprus Europe               21.8            35.1                 9.59
## 3 Croatia  Croat… Europe               24.4            36.9                 9.64
## 4 Latvia   Latvia Europe               23.6            37                  12.9 
## # ℹ 6 more variables: tobacco_code <dbl>, tobacco_levels <chr>,
## #   tobacco_custom <chr>, alcohol_code <dbl>, alcohol_levels <chr>,
## #   geom <MULTIPOLYGON [°]>

maskToClip=as.vector(st_bbox(interestingCountries))

base + geom_sf(data=interestingCountries,fill='yellow') + 
  coord_sf(xlim = c(maskToClip[1],maskToClip[3]), 
           ylim = c(maskToClip[2],maskToClip[4])) +
    geom_sf_text(data=interestingCountries,
                 aes(label=COUNTRY),
                 color='blue',
                 check_overlap = T,
                 size=3,
                 nudge_y = 0.15)

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Interactive maps

We will use leaflet library:

library(leaflet)

# create palette to fill polygons
paletteFun=colorFactor(palette = "PiYG",reverse = T, 
                       domain = countriesCIA$tobacco_levels)

#popup labels when using your cursor
popUp_labels <- sprintf("<strong>%s</strong>",
                        countriesCIA$TobaccoUse_perc) %>% lapply(htmltools::HTML)

# the base map: the WA boundaries (optional)
base_map = leaflet() %>% 
            addTiles()
# adding a layer (main map layer)
choroAll = base_map %>%
         addPolygons(data=countriesCIA,
                     stroke = F, # borders of polygon?
                     opacity =  1, # # the closer to 0 the more transparent
                     fillOpacity = 0.8, # color brigthness
                     fillColor = ~paletteFun(tobacco_levels),# coloring
                     label = popUp_labels, 
                     labelOptions = labelOptions(
                         style = list("font-weight" = "normal"),
                         textsize = "15px",
                         direction = "auto"))


# You may need a legend:
choroAll %>% addLegend(data=countriesCIA,
                        position = "bottomright",
                        pal = paletteFun,
                        values = ~tobacco_levels,
                        title = "Level",
                        opacity = 1) 

# filtering
countriesCIA_good=countriesCIA[countriesCIA$tobacco_code==0,]
countriesCIA_mean=countriesCIA[countriesCIA$tobacco_code==2,]

MyColors=c('green','orange')

# one layer per group
layer1= leaflet() %>% 
            addTiles() %>%
        addPolygons(data=countriesCIA_good,
                    color=MyColors[1],
                    fillOpacity = 1, # no transparency
                    stroke = F,
                    group = 'good') # LAYER as GROUP

layer1_2= layer1%>%addPolygons(data=countriesCIA_mean,
                               color=MyColors[2],
                               fillOpacity = 0.5, # transparency!!!
                               stroke = F,
                               group = 'average')

###### Let's add a _button_ that helps rezooming
# this is the coordinates for start location (I chose Spain)
textFun="function(btn, map){map.setView([40.416775, -3.703790], 2)}"

finalZoom= layer1_2%>%
    addEasyButton(
        easyButton(icon="fa-home", # a symbol
                   title="Zoom to Level 2", # when hovering mouse
                   onClick=JS(textFun)))

## the interactive menu
finalZoomLayers=finalZoom %>% addLayersControl(
        overlayGroups = c('good','average'),
        options = layersControlOptions(collapsed = FALSE))
finalZoomLayers

leaflet() %>% 
            addTiles() %>%
        addPolygons(data=interestingCountries) # LAYER as GROUP